!/===========================================================================/
! Copyright (c) 2007, The University of Massachusetts Dartmouth 
! Produced at the School of Marine Science & Technology 
! Marine Ecosystem Dynamics Modeling group
! All rights reserved.
!
! FVCOM has been developed by the joint UMASSD-WHOI research team. For 
! details of authorship and attribution of credit please see the FVCOM
! technical manual or contact the MEDM group.
!
! 
! This file is part of FVCOM. For details, see http://fvcom.smast.umassd.edu 
! The full copyright notice is contained in the file COPYRIGHT located in the 
! root directory of the FVCOM code. This original header must be maintained
! in all distributed versions.
!
! THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" 
! AND ANY EXPRESS OR  IMPLIED WARRANTIES, INCLUDING,  BUT NOT  LIMITED TO,
! THE IMPLIED WARRANTIES OF MERCHANTABILITY AND  FITNESS FOR A PARTICULAR
! PURPOSE ARE DISCLAIMED.  
!
!/---------------------------------------------------------------------------/
! CVS VERSION INFORMATION
! $Id$
! $Name$
! $Revision$
!/===========================================================================/

MODULE MOD_SETUP
  USE ALL_VARS
  USE MOD_PAR
  USE MOD_INPUT
  USE MOD_NCDIO
  USE MOD_NCTOOLS
  USE MOD_UTILS
  USE MOD_OBCS
  USE MOD_FORCE

  IMPLICIT NONE

  SAVE
  ! COORDINATE VARIABLES
  REAL(SP), ALLOCATABLE, TARGET :: X_GBL(:),Y_GBL(:)
  REAL(SP), ALLOCATABLE, TARGET :: X_LCL(:),Y_LCL(:)
  
  ! CORIOLIS VARIABLES
  REAL(SP), ALLOCATABLE, TARGET :: C_LCL(:)     !!CORIOLIS PARAMETER AT NODES

  ! DEPTH VARIABLE
  REAL(SP), ALLOCATABLE, TARGET :: H_LCL(:)     !! DEPTH PARAMETER AT NODES

  ! SPONGE LAYER VARIABLE
  INTEGER, ALLOCATABLE, TARGET  :: N_SPG(:)
  REAL(SP), ALLOCATABLE, TARGET :: R_SPG(:),C_SPG(:),X_SPG(:),Y_SPG(:)
  INTEGER :: NSPONGE

  PRIVATE :: SIGMA_GEOMETRIC, SIGMA_GENERALIZED, SIGMA_TANH

CONTAINS

!==============================================================================!
  SUBROUTINE SETUP_CENTER_COORDS
    USE MOD_SPHERICAL
    IMPLICIT NONE
    INTEGER I,J,IERR,STATUS,SENDID
    REAL(SP) SBUF

    INTEGER K,ITMP
    REAL(DP) VX1,VY1,VX2,VY2,VX3,VY3
    REAL(DP) EVX12,EVX13,EVX23,EVY12,EVY13,EVY23,EVXY
    REAL(DP) VX12,VY12,VX23,VY23,VX31,VY31
    INTEGER :: SENDER

    ! David Added:
    REAL(SP), allocatable :: xc_buf(:), lonc_buf(:)
    REAL(SP), allocatable :: yc_buf(:), latc_buf(:)
 
    IF(DBG_SET(DBG_SBR)) write(IPT,*) "SETUP_CENTER_COORDS: START"
    
    
!==============================================================================|
!   SET UP LOCAL MESH (HORIZONTAL COORDINATES)                                 |
!==============================================================================|
!--------------CALCULATE GLOBAL MINIMUMS AND MAXIMUMS--------------------------!
    
#  if !defined (SPHERICAL) 
! IF THIS IS NOT A SPHERICAL MODEL REMOVE MIN COORDINATE VALUE    
    
    VXMIN = MINVAL(VX(1:MT)) ; VXMAX = MAXVAL(VX(1:MT))
    VYMIN = MINVAL(VY(1:MT)) ; VYMAX = MAXVAL(VY(1:MT))
    
#  if defined (MULTIPROCESSOR) 
    IF(PAR) THEN
       !GLOBAL MIN
       SBUF = VXMIN
       CALL MPI_ALLREDUCE(SBUF,VXMIN,1,MPI_F,MPI_MIN,MPI_FVCOM_GROUP,IERR)
       SBUF = VYMIN
       CALL MPI_ALLREDUCE(SBUF,VYMIN,1,MPI_F,MPI_MIN,MPI_FVCOM_GROUP,IERR)
       !GLOBAL MAX
       SBUF = VXMAX
       CALL MPI_ALLREDUCE(SBUF,VXMAX,1,MPI_F,MPI_MAX,MPI_FVCOM_GROUP,IERR)
       SBUF = VYMAX
       CALL MPI_ALLREDUCE(SBUF,VYMAX,1,MPI_F,MPI_MAX,MPI_FVCOM_GROUP,IERR)
    END IF
#endif
    !--------------SHIFT GRID TO UPPER RIGHT CARTESIAN-----------------------------!
    
    VX = VX - VXMIN
    VY = VY - VYMIN
    VX(0) = 0.0_SP ; VY(0) = 0.0_SP
# if defined (DATA_ASSIM)
!-----------------add by Jadon Ge----------------------    
    IF(.NOT.ALLOCATED(XG))  ALLOCATE(XG(0:MGL))
    IF(.NOT.ALLOCATED(YG))  ALLOCATE(YG(0:MGL))
    IF(.NOT.ALLOCATED(XCG)) ALLOCATE(XCG(0:NGL))
    IF(.NOT.ALLOCATED(YCG)) ALLOCATE(YCG(0:NGL))
# if defined (MULTIPROCESSOR)
    IF(PAR)THEN
      SENDER = MSRID -1
      CALL MPI_BCAST(XG,MGL,MPI_F,SENDER,MPI_FVCOM_GROUP,IERR)
      CALL MPI_BCAST(YG,MGL,MPI_F,SENDER,MPI_FVCOM_GROUP,IERR)
      CALL MPI_BCAST(XCG,NGL,MPI_F,SENDER,MPI_FVCOM_GROUP,IERR)
      CALL MPI_BCAST(YCG,NGL,MPI_F,SENDER,MPI_FVCOM_GROUP,IERR)
    END IF
# endif
    XG = XG - VXMIN
    YG = YG - VYMIN
    XG(0) = 0.0_SP ; YG(0) = 0.0_SP

    XCG = XCG - VXMIN
    YCG = YCG - VYMIN
    XCG(0) = 0.0_SP ; YCG(0) = 0.0_SP
!------------------------------------------------------
# endif
    
    !--------------CALCULATE GLOBAL ELEMENT CENTER GRID COORDINATES----------------!
    CALL N2E2D(VX,XC)
    CALL N2E2D(VY,YC)
    
    XMC = XC + VXMIN
    XMC(0)= 0.0_SP
    YMC = YC + VYMIN
    YMC(0)= 0.0_SP

    IF (USE_PROJ ) THEN
        IF (SERIAL) THEN
           CALL Meters2Degrees(XMC(1:NT),YMC(1:NT)&
                &   ,PROJECTION_REFERENCE,LONC(1:NT),LATC(1:NT),NT)
# if defined (MULTIPROCESSOR)
        ELSE
            IF (MSR) THEN
                allocate(xc_buf(0:NGL),lonc_buf(0:NGL), stat=ierr)
                allocate(yc_buf(0:NGL),latc_buf(0:NGL), stat=ierr)
            END IF
            
            call ACOLLECT(MYID,MSRID,NPROCS,EMAP,XMC,xc_buf)
            call ACOLLECT(MYID,MSRID,NPROCS,EMAP,YMC,yc_buf)
            
            IF (MSR) THEN
                CALL Meters2Degrees(xc_buf(1:NGL),yc_buf(1:NGL)&
                    &   ,PROJECTION_REFERENCE,lonc_buf(1:NGL),latc_buf(1:NGL),NGL)
            END IF
            
            call ADEAL(MYID,MSRID,NPROCS,EXMAP,lonc_buf,lonc)
            call ADEAL(MYID,MSRID,NPROCS,EXMAP,latc_buf,latc)
            
            IF (MSR) THEN
                deallocate(xc_buf, lonc_buf)
                deallocate(yc_buf, latc_buf)
            END IF

# endif
        END IF

    END IF

    
# else
    ! SPHERICAL

    DO I=1,NT
       VX1=VX(NV(I,1))
       VY1=VY(NV(I,1))
       VX2=VX(NV(I,2))
       VY2=VY(NV(I,2))
       VX3=VX(NV(I,3))
       VY3=VY(NV(I,3))
       
       inner: DO K=1,1000000
          
          EVX12=VX2-VX1
          EVX13=VX3-VX1
          EVX23=VX3-VX2
          
          IF(EVX12 >  180.0_SP)THEN
             EVX12 = -360.0_SP+EVX12
          ELSE IF(EVX12 < -180.0_SP)THEN
             EVX12 =  360.0_SP+EVX12
          END IF

          IF(EVX13 >  180.0_SP)THEN
             EVX13 = -360.0_SP+EVX13
          ELSE IF(EVX13 < -180.0_SP)THEN
             EVX13 =  360.0_SP+EVX13
          END IF

          IF(EVX23 >  180.0_SP)THEN
             EVX23 = -360.0_SP+EVX23
          ELSE IF(EVX23 < -180.0_SP)THEN
             EVX23 =  360.0_SP+EVX23
          END IF
          
          EVX12=ABS(EVX12)
          EVX13=ABS(EVX13)
          EVX23=ABS(EVX23)
          
          EVY12=ABS(VY2-VY1)
          EVY13=ABS(VY3-VY1)
          EVY23=ABS(VY3-VY2)
          
          EVXY=1.E-10_SP
          
          IF((EVX12 < EVXY) .AND.(EVX13 < EVXY) .AND. (EVX23 < EVXY) &
               & .AND. (EVY12 < EVXY) .AND. (EVY13 < EVXY)  &
               & .AND. (EVY23 < EVXY)) THEN

             XC(I)=VX1
             YC(I)=VY1
             exit inner
          ELSE
             CALL ARCC(VX1,VY1,VX2,VY2,VX12,VY12)
             CALL ARCC(VX2,VY2,VX3,VY3,VX23,VY23)
             CALL ARCC(VX3,VY3,VX1,VY1,VX31,VY31)
             
             VX1=VX12
             VY1=VY12
             VX2=VX23
             VY2=VY23
             VX3=VX31
             VY3=VY31
          END IF
       END DO inner
       
    END DO

    
    XC(0) = 0.0_SP ; YC(0) = 0.0_SP

# if defined (DATA_ASSIM)
!-----------------add by Jadon Ge----------------------
# if defined (MULTIPROCESSOR)
    IF(.NOT.ALLOCATED(XG))  ALLOCATE(XG(0:MGL))
    IF(.NOT.ALLOCATED(YG))  ALLOCATE(YG(0:MGL))
    IF(.NOT.ALLOCATED(XCG)) ALLOCATE(XCG(0:NGL))
    IF(.NOT.ALLOCATED(YCG)) ALLOCATE(YCG(0:NGL))
    IF(PAR)THEN
      SENDER = MSRID -1
      CALL MPI_BCAST(XG,MGL,MPI_F,SENDER,MPI_FVCOM_GROUP,IERR)
      CALL MPI_BCAST(YG,MGL,MPI_F,SENDER,MPI_FVCOM_GROUP,IERR)
      CALL MPI_BCAST(XCG,NGL,MPI_F,SENDER,MPI_FVCOM_GROUP,IERR)
      CALL MPI_BCAST(YCG,NGL,MPI_F,SENDER,MPI_FVCOM_GROUP,IERR)
    END IF
# endif
    XG(0) = 0.0_SP ; YG(0) = 0.0_SP
    XCG(0) = 0.0_SP ; YCG(0) = 0.0_SP
!------------------------------------------------------
# endif

    LONC=XC
    LATC=YC

    IF (USE_PROJ ) THEN
        IF (SERIAL) THEN
           CALL Degrees2Meters(XC(1:NT),YC(1:NT)&
                &   ,PROJECTION_REFERENCE,XMC(1:NT),YMC(1:NT),NT)
# if defined (MULTIPROCESSOR)
        ELSE
            IF (MSR) THEN
                allocate(xc_buf(0:NGL),lonc_buf(0:NGL), stat=ierr)
                allocate(yc_buf(0:NGL),latc_buf(0:NGL), stat=ierr)
            END IF
            
            call ACOLLECT(MYID,MSRID,NPROCS,EMAP,XC,lonc_buf)
            call ACOLLECT(MYID,MSRID,NPROCS,EMAP,YC,latc_buf)
            
            IF (MSR) THEN
                CALL Degrees2Meters(lonc_buf(1:NGL),latc_buf(1:NGL)&
                    &   ,PROJECTION_REFERENCE,xc_buf(1:NGL),yc_buf(1:NGL),NGL)
            END IF
            
            call ADEAL(MYID,MSRID,NPROCS,EXMAP,xc_buf,xmc)
!JQI            call ADEAL(MYID,MSRID,NPROCS,EXMAP,xc_buf,ymc)
            call ADEAL(MYID,MSRID,NPROCS,EXMAP,yc_buf,ymc)
            
            IF (MSR) THEN
                deallocate(xc_buf, lonc_buf)
                deallocate(yc_buf, latc_buf)
            END IF

# endif
        END IF
        
    END IF


# endif

    IF(DBG_SET(DBG_SBR)) write(IPT,*) "SETUP_CENTER_COORDS: END"

  END SUBROUTINE SETUP_CENTER_COORDS
  !==============================================================================|
  !   SET HORIZONTAL MIXING_COEFFICIENT                                      |
  !==============================================================================|
  SUBROUTINE SETUP_HORIZONTAL_MIXING_COEFFICIENT
    IMPLICIT NONE
  
    if (HORIZONTAL_MIXING_KIND .eq. STTC) THEN

       if(DBG_SET(DBG_LOG)) then
          write(IPT,*) "! "
          write(IPT,*) "! Setting Staticly Variable viscosity"
       end if

       CALL LOAD_HORIZONTAL_MIXING_COEFFICIENT(NN_HVC,CC_HVC)
       
       
    else if(HORIZONTAL_MIXING_KIND .eq. CNSTNT) THEN
       
       CC_HVC=HORIZONTAL_MIXING_COEFFICIENT
       CC_HVC(0)=0.0_SP
       NN_HVC=HORIZONTAL_MIXING_COEFFICIENT
       NN_HVC(0)=0.0_SP
       
    else
       CALL FATAL_ERROR&
            &("HORIZONTAL_MIXING_KIND ERROR",&
            & "This should not happen")
       
    end if
    
    
  END SUBROUTINE SETUP_HORIZONTAL_MIXING_COEFFICIENT
  !==============================================================================|
  !   SET BOTTOM ROUGHNESS                                                       |
  !==============================================================================|
  SUBROUTINE SETUP_BOTTOM_ROUGHNESS
    IMPLICIT NONE
  
    if (BOTTOM_ROUGHNESS_KIND .eq. STTC) THEN
  
       if(DBG_SET(DBG_LOG)) then
          write(IPT,*) "! "
          write(IPT,*) "! Setting Staticly Variable Bottom Roughness"
       end if

!Changed by Siqi Li@2015/08/29
!       CALL LOAD_BOTTOM_ROUGHNESS(CC_Z0B)
       SELECT CASE(BOTTOM_ROUGHNESS_TYPE)
       CASE(BR_UDEF)
         CALL LOAD_DRAG_COEFFICIENT(CBC_UD)
       CASE DEFAULT
         CALL LOAD_BOTTOM_ROUGHNESS(CC_Z0B)
       END SELECT
!Siqi Li
              
    else if(BOTTOM_ROUGHNESS_KIND .eq. CNSTNT) THEN
       
       CC_Z0B    = BOTTOM_ROUGHNESS_LENGTHSCALE
       CC_Z0B(0) = 0.0_SP
       
    else
       CALL FATAL_ERROR&
            &("HORIZONTAL_MIXING_KIND ERROR",&
            & "This should not happen")
       
    end if
    
    
  END SUBROUTINE SETUP_BOTTOM_ROUGHNESS
  !==============================================================================|
  !   SET UP LOCAL MESH (BATHYMETRIC DEPTH)                                      |
  !==============================================================================|
  SUBROUTINE SETUP_DEPTH
    IMPLICIT NONE
    INTEGER :: IERR
    REAL(SP):: SBUF
    
    IF(DBG_SET(DBG_SBR)) write(IPT,*) "SETUP_DEPTH: START"
    HMAX = MAXVAL(H_LCL(1:MT))
    HMIN = MINVAL(H_LCL(1:MT))

# if  defined (MULTIPROCESSOR)
    IF(PAR)THEN
       SBUF = HMAX
       CALL MPI_ALLREDUCE(SBUF,HMAX,1,MPI_F,MPI_MAX,MPI_FVCOM_GROUP,IERR)
       SBUF = HMIN
       CALL MPI_ALLREDUCE(SBUF,HMIN,1,MPI_F,MPI_MIN,MPI_FVCOM_GROUP,IERR)
    END IF
# endif

    H = H_LCL

    IF(DBG_SET(DBG_SBR)) write(IPT,*) "SETUP_DEPTH: END"

  END SUBROUTINE SETUP_DEPTH

!==============================================================================|
!   SET UP LOCAL CORIOLIS FORCE                                                |
!==============================================================================|
  SUBROUTINE SETUP_CORIOLIS 
    IMPLICIT NONE
    integer:: I
    !--------------READ IN CORIOLIS PARAMETER--------------------------------------!

    IF(DBG_SET(DBG_SBR)) write(IPT,*) "SETUP_CORIOLIS: START"

    CALL N2E2D(C_LCL,COR)

#   if defined (SPHERICAL)
    IF(EQUATOR_BETA_PLANE)THEN
      F_ALFA = 1.0_SP-0.8_SP*EXP(-(COR/2.2_SP)**2)
    END IF 
#   endif

    COR = 2.*7.292e-5_SP * SIN(COR * DEG2RAD)
    
    !!  ggao for equatoral min (4deg)
!JQI20201123 -- start comment
!    IF(.NOT. EQUATOR_BETA_PLANE)THEN
!     WHERE(COR <  1.e-5_SP .AND. COR > 0.0_SP) COR =  1.e-5_SP
!     WHERE(COR > -1.e-5_SP .AND. COR < 0.0_SP) COR = -1.e-5_SP
!    END IF 
!JQI20201123 -- end comment

    IF(DBG_SET(DBG_SBR)) write(IPT,*) "SETUP_CORIOLIS: END"

  END SUBROUTINE SETUP_CORIOLIS
    
    
!==============================================================================|
!   COMPUTE GRAVITY VARIED WITH LATITUDE                                       |
!==============================================================================|

  SUBROUTINE SETUP_GRAVITY
    IMPLICIT NONE
    INTEGER I

    IF(DBG_SET(DBG_SBR)) write(IPT,*) "SETUP_GRAVITY: START"

# if defined (SPHERICAL)
    DO I=1,MT
       GRAV_N(I) = 9.780327*(1.0+0.0053024*SIN(VY(I)*DEG2RAD)**2-0.0000058*SIN(2.0*VY(I)*DEG2RAD)**2)
    END DO
    DO I=1,NT
       GRAV_E(I) = 9.780327*(1.0+0.0053024*SIN(YC(I)*DEG2RAD)**2-0.0000058*SIN(2.0*YC(I)*DEG2RAD)**2)
    END DO
# else
    GRAV_N = GRAV
    GRAV_E = GRAV
# endif        

    IF(DBG_SET(DBG_SBR)) write(IPT,*) "SETUP_GRAVITY: END"
  END SUBROUTINE SETUP_GRAVITY
    
!==============================================================================|
!   COMPUTE SPONGE LAYER FOR OPEN BOUNDARY DAMPING                             |
!==============================================================================|

  SUBROUTINE SETUP_SPONGE
    USE MOD_SPHERICAL
    IMPLICIT NONE
    REAL(SP)  TEMP,DTMP,C_SPONGE
    INTEGER :: I1, I, SENDER, IERR
    REAL(DP) X1_DP,Y1_DP,X2_DP,Y2_DP,DTMP_DP

    !--SET SPONGE PARAMETERS-------------------------------------------------------|
    
    IF(DBG_SET(DBG_SBR)) write(IPT,*) "SETUP_SPONGE: START"

    IF (NSPONGE ==0) RETURN

! NOTE: X_SPG/Y_SPG COORDINATES MUST BE AJUSTED FOR VXMIN/VYMIN

# if !defined (SPHERICAL)
    X_SPG = X_SPG - VXMIN
    Y_SPG = Y_SPG - VYMIN

    DO I=1,NT
       DO I1=1,NSPONGE
          DTMP=(XC(I)-X_SPG(I1))**2+(YC(I)-Y_SPG(I1))**2
          DTMP=SQRT(DTMP)/R_SPG(I1)

          IF(DTMP <= 1.0_SP) THEN
             C_SPONGE=C_SPG(I1)*(1.0_SP-DTMP)
             CC_SPONGE(I)=MAX(C_SPONGE,CC_SPONGE(I))
          END IF
       END DO
    END DO

# else
    ! SPHERICAL

    DO I=1,NT
       DO I1=1,NSPONGE
          X1_DP=XC(I)
          Y1_DP=YC(I)
          X2_DP=X_SPG(I1)
          Y2_DP=Y_SPG(I1)
          CALL ARC(X1_DP,Y1_DP,X2_DP,Y2_DP,DTMP_DP)
          DTMP=DTMP_DP/R_SPG(I1)

          IF(DTMP <= 1.0_SP) THEN
             C_SPONGE=C_SPG(I1)*(1.0_SP-DTMP)
             CC_SPONGE(I)=MAX(C_SPONGE,CC_SPONGE(I))
          END IF
       END DO
    END DO

# endif    

    DEALLOCATE(N_SPG,R_SPG,C_SPG,X_SPG,Y_SPG)

    IF(DBG_SET(DBG_SBR)) write(IPT,*) "SETUP_SPONGE: END"
    
    RETURN
  END SUBROUTINE SETUP_SPONGE
!==============================================================================|
  
!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%|
!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%|

!==============================================================================|
  SUBROUTINE COORDINATE_UNITS(XL,YL)
    IMPLICIT NONE
    REAL(SP), ALLOCATABLE :: XL(:),YL(:)
    integer status, ierr


    REAL(SP), allocatable :: x_buf(:), lon_buf(:)
    REAL(SP), allocatable :: y_buf(:), lat_buf(:)

    IF(DBG_SET(DBG_SBR)) write(IPT,*) "COODINATE_UNITS: START"

# if defined (SPHERICAL)

    IF (GRID_FILE_UNITS == 'degrees') THEN
       ! COPY DATA INTO VX,VY
       VX = XL
       VY = YL
       
       LON = XL
       LAT = YL
       
       IF ( USE_PROJ ) THEN
          ! USE PROJECTION TOOL BOX TO CONVERT TO METERS
        IF (SERIAL) THEN
           CALL Degrees2meters(XL(1:MT),YL(1:MT), &
               & PROJECTION_REFERENCE,XM(1:MT),YM(1:MT),MT)
# if defined (MULTIPROCESSOR)
        ELSE
            IF (MSR) THEN
                allocate(x_buf(0:MGL),lon_buf(0:MGL), stat=ierr)
                allocate(y_buf(0:MGL),lat_buf(0:MGL), stat=ierr)
            END IF
            
            call ACOLLECT(MYID,MSRID,NPROCS,NMAP,XL,lon_buf)
            call ACOLLECT(MYID,MSRID,NPROCS,NMAP,YL,lat_buf)
            
            IF (MSR) THEN
                CALL Degrees2meters(lon_buf(1:MGL),lat_buf(1:MGL)&
                    &   ,PROJECTION_REFERENCE,x_buf(1:MGL),y_buf(1:MGL),MGL)
            END IF
            
            call ADEAL(MYID,MSRID,NPROCS,NXMAP,x_buf,XM)
            call ADEAL(MYID,MSRID,NPROCS,NXMAP,y_buf,YM)
            
            IF (MSR) THEN
                deallocate(x_buf, lon_buf)
                deallocate(y_buf, lat_buf)
            END IF
          
# endif
          END IF  
        END IF
       
    ELSE IF(GRID_FILE_UNITS == 'meters') THEN

       IF ( USE_PROJ ) THEN

          ! COPY DATA INTO VX,VY
          XM = XL
          YM = YL
          
          ! USE PROJECTION TOOL BOX TO CONVERT TO METERS
        IF (SERIAL) THEN
           CALL meters2degrees(XL(1:MT),YL(1:MT), &
            & PROJECTION_REFERENCE,VX(1:MT),VY(1:MT),MT)
# if defined (MULTIPROCESSOR)
        ELSE
            IF (MSR) THEN
                allocate(x_buf(0:MGL),lon_buf(0:MGL), stat=ierr)
                allocate(y_buf(0:MGL),lat_buf(0:MGL), stat=ierr)
            END IF

            call ACOLLECT(MYID,MSRID,NPROCS,NMAP,XL,x_buf)
            call ACOLLECT(MYID,MSRID,NPROCS,NMAP,YL,y_buf)

            IF (MSR) THEN
                CALL meters2degrees(x_buf(1:MGL),y_buf(1:MGL)&
                    &   ,PROJECTION_REFERENCE,lon_buf(1:MGL),lat_buf(1:MGL),MGL)
            END IF

            call ADEAL(MYID,MSRID,NPROCS,NXMAP,lon_buf,VX)
            call ADEAL(MYID,MSRID,NPROCS,NXMAP,lat_buf,VY)

            IF (MSR) THEN
                deallocate(x_buf, lon_buf)
                deallocate(y_buf, lat_buf)
            END IF

# endif
        END IF
       
          LON = VX
          LAT = VY

       ELSE
           CALL FATAL_ERROR('You must specify a valid projection reference'&
               &,'and compile with PROJ to use files with cartesian coordinates in spherical mode')

       END IF

       
    ELSE 
       CALL FATAL_ERROR('UNRECOGNIZED GRID_FILE_UNITS: '//GRID_FILE_UNITS)
       
    END IF
    
    
# else 
!! IF NOT SPHERICAL CASE
    
    IF (GRID_FILE_UNITS == 'degrees') THEN
       
       IF (USE_PROJ) THEN
          
          LON = XL
          LAT = YL
          
          ! USE PROJECTION TOOL BOX TO CONVERT TO METERS
        IF (SERIAL) THEN
           CALL Degrees2meters(XL(1:MT),YL(1:MT), &
               & PROJECTION_REFERENCE,VX(1:MT),VY(1:MT),MT)
# if defined (MULTIPROCESSOR)
        ELSE
            IF (MSR) THEN
                allocate(x_buf(0:MGL),lon_buf(0:MGL), stat=ierr)
                allocate(y_buf(0:MGL),lat_buf(0:MGL), stat=ierr)
            END IF
            
            call ACOLLECT(MYID,MSRID,NPROCS,NMAP,XL,lon_buf)
            call ACOLLECT(MYID,MSRID,NPROCS,NMAP,YL,lat_buf)
            
            IF (MSR) THEN
                CALL Degrees2meters(lon_buf(1:MGL),lat_buf(1:MGL)&
                    &   ,PROJECTION_REFERENCE,x_buf(1:MGL),y_buf(1:MGL),MGL)
            END IF
            
            call ADEAL(MYID,MSRID,NPROCS,NXMAP,x_buf,VX)
            call ADEAL(MYID,MSRID,NPROCS,NXMAP,y_buf,VY)
            
            IF (MSR) THEN
                deallocate(x_buf, lon_buf)
                deallocate(y_buf, lat_buf)
            END IF
          
# endif
        END IF         
          XM = VX
          YM = VY
          
       ELSE
          
          CALL FATAL_ERROR('You must specify a valid projection reference'&
               &,'and compile with PROJ to use files with latitude and longitue in cartesian mode')
       END IF
       

    ELSE IF(GRID_FILE_UNITS == 'meters') THEN

       VX = XL
       VY = YL
       
       XM = XL
       YM = YL
       
       ! USE PROJECTION TOOL BOX TO CONVERT TO DEGREES
       IF (USE_PROJ ) THEN
        IF (SERIAL) THEN
           CALL meters2degrees(XL(1:MT),YL(1:MT), &
            & PROJECTION_REFERENCE,LON(1:MT),LAT(1:MT),MT)
# if defined (MULTIPROCESSOR)
        ELSE
            IF (MSR) THEN
                allocate(x_buf(0:MGL),lon_buf(0:MGL), stat=ierr)
                allocate(y_buf(0:MGL),lat_buf(0:MGL), stat=ierr)
            END IF

            call ACOLLECT(MYID,MSRID,NPROCS,NMAP,XL,x_buf)
            call ACOLLECT(MYID,MSRID,NPROCS,NMAP,YL,y_buf)

            IF (MSR) THEN
                CALL meters2degrees(x_buf(1:MGL),y_buf(1:MGL)&
                    &   ,PROJECTION_REFERENCE,lon_buf(1:MGL),lat_buf(1:MGL),MGL)
            END IF

            call ADEAL(MYID,MSRID,NPROCS,NXMAP,lon_buf,LON)
            call ADEAL(MYID,MSRID,NPROCS,NXMAP,lat_buf,LAT)

            IF (MSR) THEN
                deallocate(x_buf, lon_buf)
                deallocate(y_buf, lat_buf)
            END IF

# endif
        END IF
          
       END IF

    ELSE 
       CALL FATAL_ERROR('UNRECOGNIZED GRID_FILE_UNITS: '//TRIM(GRID_FILE_UNITS))
       
    END IF

# endif

    IF(DBG_SET(DBG_SBR)) write(IPT,*) "COODINATE_UNITS: END"
    
  END SUBROUTINE COORDINATE_UNITS

!==============================================================================|
! SETUP THE SIGMA COORDINATES FOR THE MODEL                                    |
!==============================================================================|
!==============================================================================|
! This program is used to set up the coordinate in the vertical.               !
!								               !
! case(1): sigma levels                                                        !
! sigma levels are determined by a formula of                                  !
!                      sigma(k)=-[(k-1)/(kb-1)]^k11                            !
!    p_sigma=1: uniform sigma layers                                           !
!    p_sigma=2: layers satisfying a parabolic function with high               !
!               vertical resolution near the surface and bottom.               !
!    p_sigma can be used any real number                                       !
!									       !
! case(2): general vertical level                                              !
! vertical levels are determined by the formula                                !
!                tanh[(dl+du)((kbm1-k)/kbm1)-dl]+tanh(dl)                      !
!        z(k)= ------------------------------------------  - 1                 !
!                      tanh(dl) + tanh(du)                                     !
!                                                                              !
! case(3): constant layer transformation                                       !
! four values need to be specified:                                            !
!  DUU the upper boundaries, up to which the co-ordinates are parallel must be !
!      defined.                                                                !
!  DLL the lower boundaries, up to which the co-ordinates are parallel must be !
!      defined.                                                                !
!  HMIN1 the minimum water depth at which the layers are constant. If H < HMIN1!
!      then sigma co-ordinates are used.                                       !
!                                                                              !
! Reference of case(2), case(3) and case(4):                                   !
! Pietrzak, J.D., Jan B. Jakobson, Hans Burchard, Hans Jacob Vested, Ole       !
! Petersen , 2002. A three-dimensional hydrostatic model for coastal and ocean !
! modelling using a generalised topography following co-ordinate system. Ocean !
! Modelling 4, 173-205                                                         !
!                                                                              !
!  calculates: z(m,kb) vertical levels					       !
!  calculates: dz(m,kb-1) delta between vertical levels		               !
!  calculates: zz(m,kb-1) intra-vertical levels				       !
!  calculates: dzz(m,kb-2) delta between intra-vertical levels		       !
!==============================================================================|

  SUBROUTINE SETUP_SIGMA          
    !==============================================================================|
    IMPLICIT NONE
    INTEGER :: K,KK
    INTEGER :: I
    REAL(SP):: ZTMP(KB)
    REAL(SP):: X1,X2,X3  
    REAL(SP):: DR,RCL,RCU
    !==============================================================================|

    IF(DBG_SET(DBG_SBR)) WRITE(IPT,*)"SETUP_SIGMA: START"

    IF(DBG_SET(DBG_SBRIO)) THEN
       WRITE(IPT,*)"==================="
       WRITE(IPT,*)"  SET_SIGMA IO"
       WRITE(IPT,*)"  STYPE = "//TRIM(STYPE)
       WRITE(IPT,*)"  P_SIGMA = ", P_SIGMA
       WRITE(IPT,*)"  KB = ", KB
       WRITE(IPT,*)"==================="
    END IF

    !---------------------------------------------
    !---------------------------------------------
    SELECT CASE(STYPE)
    !---------------------------------------------
    !SIGMA_COORDINATE_TYPE = UNIFORM (DEGENERATE CASE OF GEOMETRIC)
    CASE(STYPE_UNIFORM)
    !---------------------------------------------
       IF(P_SIGMA > 1 .AND. MOD(KB,2) == 0) &
            CALL FATAL_ERROR('SETUP_SIGMA: COORDINATE TYPE:'//trim(STYPE)&
            &,'kb shoude be an odd number for this type of sigma coordinates....' )
       CALL SIGMA_GEOMETRIC
    !---------------------------------------------
    !SIGMA_COORDINATE_TYPE = GEOMETRIC
    CASE(STYPE_GEOMETRIC)
    !---------------------------------------------

       IF(P_SIGMA > 1 .AND. MOD(KB,2) == 0) &
            CALL FATAL_ERROR('SETUP_SIGMA: COORDINATE TYPE:'//trim(STYPE)&
            &,'kb shoude be an odd number for this type of sigma coordinates....' )
       CALL SIGMA_GEOMETRIC
    !---------------------------------------------
    ! SIGMA_COORDINATE_TYPE = TANH
    CASE(STYPE_TANH)
    !---------------------------------------------   
       CALL SIGMA_TANH
    !---------------------------------------------   
    !SIGMA_COORDINATE_TYPE = GENERALIZED
    CASE(STYPE_GENERALIZED)
    !---------------------------------------------   
       CALL SIGMA_GENERALIZED

    ! THIS IS A CURRENTLY UNUSED METHOD
!!$!---------------------------------------------   
!!$!SIGMA_COORDINATE_TYPE = WHAT THE HELL IS THIS?
!!$   CASE("UNKNOWN")
!!$!---------------------------------------------   
!!$
!!$    CALL SIGMA_UNKNOWN
    CASE DEFAULT
       CALL FATAL_ERROR("SET_SIGMA: REACHED DEFAULT CASE FOR SIGMA COOR&
            &DINATE TYPE")
    END SELECT

    !---------COMPUTE SIGMA DERIVATIVES AND INTRA SIGMA LEVELS---------------------!
#    if defined (MULTIPROCESSOR)
    IF(PAR)CALL AEXCHANGE(EC,MYID,NPROCS,Z1)
    IF(PAR)CALL AEXCHANGE(NC,MYID,NPROCS,Z)
#    endif

    IF(DBG_SET(DBG_SBR)) WRITE(IPT,*)"SETUP_SIGMA: END"

    RETURN

  END SUBROUTINE SETUP_SIGMA
!==============================================================================|
  SUBROUTINE SIGMA_GEOMETRIC
      IMPLICIT NONE
      INTEGER :: I,K
      REAL(SP):: ZTMP(KB)
      !orginal formula to set sigma   
!      IF(P_SIGMA == 1)THEN
      IF(ABS(P_SIGMA - 1.0) < 1.0E-6)THEN
         DO K=1,KB
!JQI NOV2021
!            ZTMP(K) = -((K-1)/FLOAT(KB-1))**P_SIGMA 
            ZTMP(K) = -((K-1)/FLOAT(KB-1))**P_SIGMA_DISTRIBUTION        !PNNL added this option P_SIGMA_DISTRIBUTION 
!JQI NOV2021
        END DO
      ELSE
         DO K=1,(KB+1)/2
            ZTMP(K) = -((K-1)/FLOAT((KB+1)/2-1))**P_SIGMA/2 
         END DO
         DO K=(KB+1)/2+1,KB
            ZTMP(K) = ((KB-K)/FLOAT((KB+1)/2-1))**P_SIGMA/2-1.0
         END DO
      END IF

      DO I=1,M
         DO K=1,KB
            Z(I,K)=ZTMP(K)
         END DO
      END DO

      DO I=1,N
         DO K=1,KB
            Z1(I,K)=(Z(NV(I,1),K)+Z(NV(I,2),K)+Z(NV(I,3),K))/3.0_SP
         END DO
      END DO
    END SUBROUTINE SIGMA_GEOMETRIC
!--------------------------------------------------------------------    
    SUBROUTINE SIGMA_GENERALIZED
      IMPLICIT NONE
      INTEGER :: I,K, kk
      REAL(SP):: X1,X2,X3  
      REAL(SP):: DR,RCL,RCU
 
      DO I=1,M
         IF(H(I) < HMIN1)THEN
            Z(I,1)=0.0
            DL2=0.001;DU2=0.001
            DO K=1,KBM1
               X1=DL2+DU2
               X1=X1*(KBM1-K)/KBM1
               X1=X1-DL2
               X1=TANH(X1)
               X2=TANH(DL2)
               X3=X2+TANH(DU2)

               Z(I,K+1)=(X1+X2)/X3-1.0_SP
            END DO
         ELSE
            DR=(H(I)-DUU-DLL)/H(I)/(KB-KU-KL-1)

            Z(I,1)=0.0_SP

            DO K=2,KU+1
               Z(I,K)=Z(I,K-1)-ZKU(K-1)/H(I)
            END DO

            DO K=KU+2,KB-KL
               Z(I,K)=Z(I,K-1)-DR
            END DO

            KK=0
            DO K=KB-KL+1,KB
               KK=KK+1
               Z(I,K)=Z(I,K-1)-ZKL(KK)/H(I)
            END DO
         END IF
      END DO

      DO I=1,N
         DO K=1,KB
            Z1(I,K)=(Z(NV(I,1),K)+Z(NV(I,2),K)+Z(NV(I,3),K))/3.0_SP
         END DO
      END DO
    END SUBROUTINE SIGMA_GENERALIZED
!--------------------------------------------------------------------
    SUBROUTINE SIGMA_TANH
      IMPLICIT NONE
      INTEGER :: I,K
      REAL(SP):: X1,X2,X3  
    
      Z=0.0;Z1=0.0
      DO K=1,KBM1
         X1=DL2+DU2
         X1=X1*(KBM1-K)/KBM1
         X1=X1-DL2
         X1=TANH(X1)
         X2=TANH(DL2)
         X3=X2+TANH(DU2)
         DO I=1,M
            Z(I,K+1)=(X1+X2)/X3-1.0_SP
         END DO
         DO I=1,N
            Z1(I,K+1)=(X1+X2)/X3-1.0_SP
         END DO
      END DO
    END SUBROUTINE SIGMA_TANH
!--------------------------------------------------------------------
!!$ SUBROUTINE SIGMA_UNKNOWN
!!$   DO I=1,M
!!$     IF(H(I) < HMIN1)THEN
!!$       RCU=-DUU/HMIN1
!!$       RCL=DLL/HMIN1-1
!!$       DR=(RCL-RCU)/(KB-KU-KL-1)
!!$
!!$       DO K=1,KU
!!$         ZKU(K)=RCU/KU
!!$       END DO
!!$       DO K=1,KL
!!$         ZKL(K)=(-1.0_SP-RCL)/KL
!!$       END DO
!!$
!!$       Z(I,1)=0.0_SP
!!$
!!$       DO K=2,KU+1
!!$         Z(I,K)=Z(I,K-1)+ZKU(K-1)
!!$       END DO
!!$
!!$       DO K=KU+2,KB-KL
!!$         Z(I,K)=Z(I,K-1)+DR
!!$       END DO
!!$
!!$       KK=0
!!$       DO K=KB-KL+1,KB
!!$         KK=KK+1
!!$         Z(I,K)=Z(I,K-1)+ZKL(KK)
!!$       END DO
!!$
!!$     ELSE
!!$       DR=(H(I)-DUU-DLL)/H(I)/(KB-KU-KL-1)
!!$
!!$       Z(I,1)=0.0_SP
!!$
!!$       DO K=2,KU+1
!!$         Z(I,K)=Z(I,K-1)-ZKU(K-1)/H(I)
!!$       END DO
!!$
!!$       DO K=KU+2,KB-KL
!!$         Z(I,K)=Z(I,K-1)-DR
!!$       END DO
!!$
!!$       KK=0
!!$       DO K=KB-KL+1,KB
!!$         KK=KK+1
!!$         Z(I,K)=Z(I,K-1)-ZKL(KK)/H(I)
!!$       END DO
!!$     END IF
!!$   END DO
!!$   
!!$   DO I=1,N
!!$     DO K=1,KB
!!$       Z1(I,K)=(Z(NV(I,1),K)+Z(NV(I,2),K)+Z(NV(I,3),K))/3.0_SP
!!$     END DO
!!$   END DO    
!!$END SUBROUTINE SIGMA_UNKNOWN

!==============================================================================|

  SUBROUTINE SETUP_SIGMA_DERIVATIVES
    USE ALL_VARS
    IMPLICIT NONE
    INTEGER :: K, I


    IF(DBG_SET(DBG_SBR)) WRITE(IPT,*)"SETUP_SIGMA_DERIVATIVES: START"

    DO K=1,KB-1
       DO I=1,MT
          DZ(I,K)  = Z(I,K)-Z(I,K+1)
          ZZ(I,K)  = .5_SP*(Z(I,K)+Z(I,K+1))
       END DO
       DO I=1,NT
          DZ1(I,K)  = Z1(I,K)-Z1(I,K+1)
          ZZ1(I,K)  = .5_SP*(Z1(I,K)+Z1(I,K+1))
       END DO
    END DO

    DO I=1,MT
       ZZ(I,KB) = 2.0_SP*ZZ(I,KB-1)-ZZ(I,KB-2)
    END DO
    DO I=1,NT
       ZZ1(I,KB) = 2.0_SP*ZZ1(I,KB-1)-ZZ1(I,KB-2)
    END DO

    DO K=1,KBM2
       DO I=1,MT
          DZZ(I,K) = ZZ(I,K)-ZZ(I,K+1)
       END DO
       DO I=1,NT
          DZZ1(I,K) = ZZ1(I,K)-ZZ1(I,K+1)
       END DO
    END DO

    DZZ(:,KBM1) = 0.0_SP
    DZ(:,KB)    = 0.0_SP
    DZZ1(:,KBM1) = 0.0_SP
    DZ1(:,KB)    = 0.0_SP


    !----------OUTPUT VALUES-TO INFOFILE-------------------------------------------!

    IF(DBG_SET(DBG_LOG)) THEN
       WRITE(IPT,*  )'!'
       WRITE(IPT,*  )'!'
       WRITE(IPT,*)'!                SIGMA LAYER INFO     '
       WRITE(IPT,*) "SIGMA TYPE:",TRIM(STYPE)
       WRITE(IPT,70)
       SELECT CASE(STYPE)
       CASE(STYPE_UNIFORM)
          DO K=1,KB
             WRITE(IPT,80) K,Z(1,K),ZZ(1,K),DZ(1,K),DZZ(1,K)
          END DO
       CASE(STYPE_RESTART)
          DO K=1,KB
             WRITE(IPT,80) K,Z(1,K),ZZ(1,K),DZ(1,K),DZZ(1,K)
          END DO
       CASE(STYPE_GEOMETRIC)
          DO K=1,KB
             WRITE(IPT,80) K,Z(1,K),ZZ(1,K),DZ(1,K),DZZ(1,K)
          END DO
       CASE(STYPE_TANH)
          DO K=1,KB
             WRITE(IPT,80) K,Z(1,K),ZZ(1,K),DZ(1,K),DZZ(1,K)
          END DO
       CASE(STYPE_GENERALIZED) ! THIS IS CASE SPECIFIC 
          WRITE(IPT,*)"SET CASE SPECIFIC GENERALIZED SIGMA LAYER OUTPUT TO SCREEN &
               &IN mod_setup.F"
       END SELECT
       WRITE(IPT,*  )'!'
    END IF

    IF(DBG_SET(DBG_SBR)) WRITE(IPT,*)"END SETUP_SIGMA_DERIVATIVES"

    !----------FORMAT STATEMENTS---------------------------------------------------!

70  FORMAT(2x,'k',13x,'z',11x,'zz',11x,'dz',11x,'dzz')
80  FORMAT(' ',i5,4f13.8)


    IF(DBG_SET(DBG_SBR)) WRITE(IPT,*)"SETUP_SIGMA_DERIVATIVES: END"
  END SUBROUTINE SETUP_SIGMA_DERIVATIVES
  
END MODULE MOD_SETUP

